keep<-grep("_share",colnames(ag))
keep<-c(keep, which(colnames(ag)=="pop"), which(colnames(ag)=="state"), which(colnames(ag)=="yes"),which(colnames(ag)=="treat"), which(colnames(ag)=="median.age."), which(colnames(ag)=="median.age..male.population"),which(colnames(ag)=="median.household.income..in.2017.inflation.adjusted.dollars."), which(colnames(ag)=="population.density..per.sq..mile."), which(colnames(ag)=="two_party_trump_pct") )


ag2<-ag[,keep]
ag2$weight<-ag2$pop

start<-1
stop<-which(colnames(ag2)=="tot_clr_index_property5yr_share")


est_vars<-colnames(ag2)[start:stop]
est_vars<-c("pop",est_vars, "median.age." ,"median.age..male.population" ,"median.household.income..in.2017.inflation.adjusted.dollars.","two_party_trump_pct" )
est_vars

d<-ag2

split <- function(value){
	median <- median(value, na.rm=T)
	for(j in 1:length(value)){
		if(!is.na(value[j])){
			if(value[j] <= median){
				value[j] <- "Low"
			}else{
				value[j] <- "High"
			}
		}
	}
	return(value)
}
	
for(i in 1:length(est_vars)){
	d[,est_vars[i]] <- split(d[,est_vars[i]])
	
	
	#d[,est_vars[i]]<-(d[,est_vars[i]] - mean(d[,est_vars[i]], na.rm=T)) / (sd(d[,est_vars[i]], na.rm=T))
	
}

coef<-NA
se<-NA
lb<-NA
ub<-NA
sig<-NA
p.val<-NA
treat<-NA

d2<- d[d$treat==1,]

for(i in 1:length(est_vars)){
	m<-lm(yes~ d2[,est_vars[i]], data=d2, weights=d2$weight)
	coef[i]<-m$coefficients[2]
	
	m2<-coeftest(m, vcov=vcovHC(m))
	
	
	se[i]<-m2[2,"Std. Error"]
	lb[i]<-coef[i]-1.96*se[i]
	ub[i]<-coef[i]+1.96*se[i]
	sig[i]<-I(lb[i]*ub[i]>0)
	p.val[i]<-m2[2,"Pr(>|t|)"]
	treat[i] <- "Control"
	
}

d3 <- d[d$treat==2 & d$state != "NJ",]
for(i in 1:length(est_vars)){
	m<-lm(yes~ d3[,est_vars[i]], data=d3, weights=d3$weight)
	coef[i+44]<-m$coefficients[2]
	
	m2<-coeftest(m, vcov=vcovHC(m))
	
	se[i+44]<-m2[2,"Std. Error"]
	lb[i+44]<-coef[i]-1.96*se[i]
	ub[i+44]<-coef[i]+1.96*se[i]
	sig[i+44]<-I(lb[i]*ub[i]>0)
	p.val[i+44]<-m2[2,"Pr(>|t|)"]
	treat[i+44] <- "Treatment"
	
}

res<-cbind.data.frame(var=est_vars,coef, se, lb, ub,p.val,sig, treat)

res<-res[order(res$treat, res$sig,decreasing=F),]

p <- p.adjust(p=res$p.val, method = "BH", n = length(res$p.val))
res$p2<-p
res$sig2<-I(res$p2<.05)
res<-res[order(res$sig2,decreasing=T),]

res$sig3 <- "P>.05"
res$sig3[res$sig2==T] <- "P<.05"

res2 <- read.csv("data/correlation_names.csv")
res3 <- left_join(res,res2)

res3$coef <- res3$coef*100
res3$var[is.na(res3$label)]
res3 <- na.omit(res3)
cor <- res3 %>%
	arrange(treat,coef) %>%    # First sort by val. This sort the dataframe but NOT the factor levels
	mutate(label=factor(label, levels=unique(label)))   # This trick updates the factor levels

ggplot(aes(y=coef,x=label, color=sig3), data = cor[cor$treat== "Control",])+
	geom_point(position = position_dodge(width=.5))+
	ylim(-15,15)+
	theme_bw()+
	xlab("")+
	ylab("Percentage-point change in agreement from a one\n standard deviation increase in the predictor") + 
	geom_hline(yintercept=0, linetype="dotted") + 
	theme(panel.spacing = unit(1, "lines")) + 
	theme(panel.grid.major = element_line(colour = "white"),panel.grid.minor = element_line(colour = "white"), axis.title.x = element_text(vjust=-0.5)) +
	theme(legend.position="bottom") +
	labs(col = "Benjamini-Hochberg corrected p-value") + 
	#theme(axis.text.x = element_text(angle=90,hjust=0.95)) +
	theme(axis.title.y = element_text(hjust=1)) +
	#ggtitle("Response to collaboration discussion requests") +
	coord_flip() #+
	#facet_grid(.~treat)
cor
ggsave("results/figure2.pdf",width=8, height=7.5)